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ABSTRACT 

We investigate the dynamical evolution of a star cluster in an external tidal field by using N- 
body simulations, with focus on the effects of the presence or absence of neutron star natal 
velocity kicks. We show that, even if neutron stars typically represent less than 2% of the total 
bound mass of a star cluster, their primordial kinematic properties may affect the lifetime of 
the system by up to almost a factor of four. We interpret this result in the light of two known 
modes of star cluster dissolution, dominated by either early stellar evolution mass loss or two- 
body relaxation. The competition between these effects shapes the mass loss profile of star 
clusters, which may either dissolve abruptly (“jumping”), in the pre-core-collapse phase, or 
gradually (“skiing”), after having reached core collapse. 
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1 INTRODUCTION 


The rate at which star clusters lose mass has been one of the en- 
during problems of st ellar dynamics. In one of the earliest results, 
lAmbartsumianI l ll938h already highlighted the role of relaxation. 
The landmark survey of ichernoff & Weinberd ( Il99(lh added a mass 
spectrum, stellar evolution and a tidal boundary, and also revealed 
the importance of the initial structure of the star cluster. But several 
other factors also influenc e the lifetime of star clusters, including 
the binary population (e.g.lTanikawa & Fukushi gell2009ti . the form 
of the Galactic orbit fe.g. lBaumgardt & Makinc 2003|), the form of 
the G alactic potential and tidal sh ocking (e.g. Gnedin & Ostriken 
Il997h . and the crossing time scale l lWhitehead et alj|2013h . 


In this Letter we add one more influence: natal kicks of neutron 
stars (NS). Though neutron stars may account for less than 2% of 
the cluster by mass, we find, astonishingly, that the presence or 
absence of kicks may change the lifetime of a star cluster by almost 
a factor of four. Though the existence of natal kicks of neutron 
stars is not in doubt, their d istribution and dispersion a re difficult to 
establish (see, for example. IPodsiadlowski et ^11200511 . 


In order to isolate the effect of this one factor we consider mod¬ 
els fro m another landmark survey of the evolution of star clusters: 
that bv iBaumgardt & Makind ( 120031 hereafter BM03). As it hap¬ 
pens, they imposed no natal kicks on neutron stars, and it was the 
attempt to reproduce some of their results that led to our discov¬ 
ery. Indeed their principal models, which begin with a King profile 
with Wo = 5, evolve very differently, both qualitatively and quan¬ 
titatively, if natal kicks are applied. 

The particular models we considered are described in the fol¬ 
lowing section, while Sect. presents our results in some detail, 
including some information on core collapse and mass segregation. 


The final section summarises our conclusions, and attempts to in¬ 
terpret them in the context of other recent work. 


2 DESCRIPTION OF THE RUNS 


We simulate the evolution of a globular cl uster as in IBM03L but 
using NBODY6 jNitadori & Aarsetlj|2012ll . We have performed a 
survey of simulations in an accelerating, non-rotating frame, using 
a number of particles between N = 8192 and N = 131072, a 
Kroupa IMF dKroupalboOlh . with the mass of the stars between 
0.1 Mq and 15 Mq (resulting in a theoretical mean mass (m) = 
0.547 Mq), and metallicity Z = 0.001. Natal kicks, when they 
were applie d, had a Maxwellian distri bution with a = 190 km s“^ 
(see eq. 3 in iHansen & Phinn^ll997h . 

In our simulations, the cluster is in a circular orbit, or in an el¬ 
liptical orbit with eccentricity e = 0.5, in a logarithmic Galactic 
potential 0 = Vg^ ln(i?G), where Vg is the circular velocity and 
Rg is the Galactocentric dista nce. For the majority of our runs we 
have used a Roche-lobe filling lKin 3 (Il96d) model with Wo = 5 as 
initial condition. Additional simulations have been performed by 
increasing the initial concentration of the King profile {Wo = 7). 
The clusters start at a Galactic radius of 8.5 kpc, with an initial ve¬ 
locity of 220 km s“^ (in the circular case); in the elliptical case the 
apogalacticon is at 8.5 kpc and the initial speed there is reduced ap¬ 
propriately, while the size of the cluster is determined by assuming 
a Roche-lobe filling condition at perigalacticon. The initial condi- 
tions for all the sim ulations have been generated using McLuster 
dKiipper et aklbOlih . 

The tidal radius of the cluster was defined as the Jacobi radius 
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Table 1. N -body simulation properties 


Model 

N 

Wo 

e 

Mo 

[M©] 

Th 

[pc] 

rj 

[pc] 

'^diss 

[Myr] 

rpBM 

diss 

[Myr] 

Tec 

[Myr] 

rpBM 

[Myr] 

8kK 

8192 

5.0 

0.0 

4497.3 

4.53 

24.35 

2426 

- 

2666 

- 

16kK 

16384 

5.0 

0.0 

8990.7 

5.73 

30.67 

2816 

- 

- 

- 

32kK 

32768 

5.0 

0.0 

18419.2 

7.23 

38.96 

3669 

- 

- 

- 

64kK 

65536 

5.0 

0.0 

36183.1 

9.10 

48.79 

4516 

- 

- 

- 

128kK* 

131072 

5.0 

0.0 

71422.0 

11.46 

61.21 

5927 

- 

- 

- 

8kN 

8192 

5.0 

0.0 

4497.3 

4.53 

24.35 

4137 

4149 

3142 

3329 

16kN 

16384 

5.0 

0.0 

8990.7 

5.73 

30.67 

5932 

6348 

4810 

5062 

32kN 

32768 

5.0 

0.0 

18419.2 

7.23 

38.96 

9384 

9696 

7788 

8412 

64kN 

65536 

5.0 

0.0 

36183.1 

9.10 

48.79 

14414 

15197 

12375 

13193 

128kN 

131072 

5.0 

0.0 

71659.0 

11.46 

61.27 

22707 

23769 

20307 

21339 

128kKe 

131072 

5.0 

0.5 

71453.0 

5.50 

29.43 

5479 

- 

6859 

- 

128kNe 

131072 

5.0 

0.5 

71453.0 

5.50 

29.43 

11254 

11675 

8952 

9332 

128kK7 

131072 

7.0 

0.0 

71780.9 

7.14 

61.31 

18369 

- 

18267 

- 

128kN7 

131072 

7.0 

0.0 

71780.9 

7.14 

61.31 

24494 

25506 

11886 

12620 


Note. — The capital letter in the model label indicates if the model is characterized by the presence (#K, e.g. 
128kK) or the absence (#N, e.g. 128kN) of NS initial kicks. The star (*) denotes a model for which two different 
numerical realizations have been evolved; the values are the average of those for the two simulations. 


where M is the “hound” cluster mass. The quantities M and rj 
were determined self-consistently and iteratively by first assuming 
that all stars are still bound and calculating the tidal radius with this 
formula. In a second step, we calculated the mass of all stars inside 
r j relative to the density centre of all stars, and used it to obtain a 
new estimate for r j. This method was repeated until convergence. 
Escapers were not removed from the simulations. 

The properties of the simulations are presented in Table [T] The 
significance of the model label is stated in the note to the Table. 
Column 4 gives the orbital eccentricity; columns 5, 6 and 7 are the 
initial values of the total bound mass, the half-mass radius and the 
Jacobi radius, re spective ly. Column 8 gives the dissolution time, 
which, following IbMOJI . is defined as the time when 95% of the 
mass was lost from the cluster, while column 10 gives th e core¬ 
collapse time. The corresponding quantities from IbMOJI are re¬ 
ported in columns 9 and 11, respectively. In our analysis, the mo¬ 
ment of core collapse Tec has been determined by inspecting the 
time evolution of the core radius and of the innermost lagrangian 
radius enclosing 1% of the total mass. 


3 RESULTS 

3.1 Lifetime and main properties of the models 

The main result of our investigation is that the presence or the ab¬ 
sence of NS natal velocity kicks can affect significantly the lifetime 
of star clusters, up to almost a factor of four. This striking result is 
illustrated by a series of “reference models” (with Wo = 5, e = 0, 
and N = 128k, 64k, 32k, 16k, and 8k), with or without NS velocity 
kicks; the time evolution of the bound mass of these models is pre¬ 
sented in Fig.[T] The difference in the behaviour of the models with 
or without NS kicks starts early in their evolution {M/Mo ~ 0.8) 
and leads to a dramatic contrast in the slope of the graph at the final 
stages of evolution {M/Mq < 0.2). 

An important aspect of the very rapid dissolution of the models 
with NS kicks is that, in almost all cases, they fail to reach core col¬ 
lapse during their evolution, as opposed to the models without NS 


kicks, which show signatures of core collapse at a time correspond¬ 
ing to 0.1 < M/Mo < 0.2. The only exception is given by model 
8kK, which reaches core collapse at the very late stages of evolu¬ 
tion, 240 Myr after the formal dissolution time (see Tab.|T] row 1, 
and the corresponding black dot in Fig.|T]l. The fact that it reaches 
core collapse at all, while larger models do not, is attributable to its 
short initial relaxation time. 

We have also considered two additional pairs of models, as 
representative cases of the regime of high initial concentration 
(Wo = 7; models 128kK7 and 128kN7) and of the evolution of 
a star cluster on an elliptic orbit (e = 0.5; models 128kKe and 
128kNe). Here the effects of the presence of NS kicks on the star 
cluster lifetime are less severe compared to those on the “refer¬ 
ence models”, but they are still significant (see Fig.j^. Both mod¬ 
els 128kK7 and 128kKe reach core collapse, although at a very 
late stage of evolution. Of the systems without NS kicks, model 
128kNe reaches core collapse at a mass comparable to that of “ref¬ 
erence models” without NS kicks, while model 128kN7 has the 
largest mass at Tec of all the models in our survey which reach core 
collapse; such a result is not surprising, given its initial concentra¬ 
tion. 

Another useful diagnostic of the differences between models 
with and without NS kicks is provided by the mean mass of stars 
in the innermost lagrangian shell, enclosing 1% of the total bound 
mass of a system. Its time evolution is illustrated in Fig. 0 for all 
models in our survey with N = 128k particles. In almost all cases, 
the mean mass in the innermost shell initially shows a decrease, 
which is due to the early evolution and escape of massive stars; as 
expected, this effect is more pronounced for systems with NS kicks. 
Nonetheless, after only a few Gyr, the value of the central mean 
mass starts to increase, reflecting the process of mass segregation. 
For models that reach core collapse, the mean mass in the final 
stages of evolution falls within the range 1.2 ^ (m) ^ 1.4Mq, 
which indicates the dominance of neutron stars in the central re¬ 
gions of the system. Not surprisingly, the rapidly dissolving model 
128kK (Fig.[^ red line) shows a final mean mass which is compa¬ 
rable to the initial value. 
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Figure 1. Time evolution of the fraction of bound mass (normalized to the 
initial value) of models with initial concentration Wq = 5, on circular 
orbits. The models are characterized by different number of particles and 
by the presence (red lines; right to left: 128kK, 64kK, 32kK, 16kK and 
8kK) or the absence (blue lines; right to left: 128kN, 64kN, 32kN, 16kN 
and 8kN) of NS kicks. The black d ots mar k when core collapse occurs. The 
con'esponding models studied bv IbMOsI i.e. without NS kicks, are also 
shown (green dashed lines; right to left: 128k, 64k, 32k, 16k and 8k); this 
data was retrieved by means of the data extraction tool Dexter. 


3.2 Detailed comparison with Baumgardt & Makino (2003) 


Despite our best efforts in repro ducing the initial conditions and the 
numerical set-up described bv iBMQSl we note that there are still 
non-negligible discrepancies between our models without NS natal 
kicks and the corresponding ones in their original investigation (see 
Table[T]and Figs. [T] and [^. We have attempted to identify the main 
reasons for these discrepancies in the intrinsic differences between 
the W-body codes used to perform the simulations, and in particular 
slightly different stellar evolution prescriptions. 

We per formed our simulations b y usin g the G PU version of 
NBODY6 dNitadori & Aarsethl2012h. whil e BMO^ used the public 
GRAPE-6 version of NBODY4 lAarsethll 19991) . The latter treats 
components of binaries as single stars, without collisions or ex¬ 
change of mass, and the resulting differences might partially ex¬ 
plain the increasing discrepancy after core collapse for the models 
depicted in Fig.[2 bec ause of t he increase in the number of binaries 
at this time. Moreover, BMQ3| used a presc ription for the properties 
of stellar remnants by Hurley et^ ( l2000h . while in NBODY6 the 
lEldridge & ToutI 1 I2OO4I) recipe is now use d. To test this , we ca rried 
out a simulation of model 128kN with the iHurlev et alj 1 I 2 OOOI) pre¬ 
scription for stellar remnants, but we obtained a dissolution time of 
Tdiss = 23.0 Gyr, which reduces the discrepancy by only about 
30% (see data for model 128kN in Table[TJ. 

To assess stochastic effects (such as run-to-run variations) we 
also performed additional simulations of models 128kN and 64kN 
by evolving different numerical realizations of the same initial con¬ 
ditions, and by evolvin g the sa me realization in several indepen¬ 
dent simulations (as in IbMOSH . Finally, we performed a simula¬ 
tion of model 128kN in which the escapers were progressively re¬ 
moved (as in iBMOSi) . but again without any significant difference 
(Tdias = 22.9 Gyr). 


Figure 2. Time evolution of the fraction of bound mass of models with 
(i) Wq = 5, on elliptic orbits; (ii) Wq = 7, on circular orbits. As in 
Fig.O models with NS kicks are denoted by red lines (right to left: 128kK7 
and 128kKe), and without NS kicks by blue lines (right to left: 128kN7 
and 128kNe). Dashed green lines show the corresponding models (without 
kicks) from BM03. 



t [Myr\ 


Figure 3. Evolution of the mean mass of the stars in the innermost lan- 
grangian shell, containing 1 % of the bound mass, evaluated for all models 
with N=128k. The vertical arrows mark the moment of core collapse (in the 
five models which exhibit core collapse). 


None of these effects was able, individually, to account for 
the observed discrepancy. Therefore, we believe that the small but 
systematic discrepancy be tween o ur models without NS kicks and 
the corresponding ones in iBMOSi results from a combination of all 
the effects mentioned above, and others which we have not stud¬ 
ied, including possible differences in the way in which models are 
virialised and scaled in different codes. As we shall show later 
(Sect. l4Tt the sensitivity of these runs to small effects is such that 
apparently trivial differences could have significant effects. 
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4 DISCUSSION 

4.1 Two modes of star cluster dissolution 

We have found that the presence or absence of neutron star kicks, in 
the models we have studied, can change the lifetime of a star cluster 
by a large factor. We shall now try to interpret our results in the 
context of previous studies of star cluster dissolution mechanisms, 
with the aim of understanding why it is that a process which affects 
such a small fraction of the mass can have such a dramatic effect. 

We consider initially tidally filling, multi-mass models with 
stellar evolution. Over the years, several numerical investiga¬ 
tions have shown that the dissolution time is strongly affected 
by two factors: the initial relaxation time and the initial con- 
centration (represented by th e King parameter Wo)- In particular 
IChernoff & Weinberd l ll990l) showed that, for a Salpeter-like IMF, 
their models with Wo = 1 or 3 dissolved quickly, in less than a Gyr, 
and without core collapse, while models with Wo = 7 all entered 
core collapse, after about 10 Gyr or longer. Clusters with Wo = 3 
and a steeper IMF (and hence a longer time scale for mass loss 
by stellar evolution) could enter core collapse before dissolution, 
provided that relaxation was fast enough. Thus there is a tension 
between the time scales of stellar evolution and relaxation, which 
plays out di fferently depending on th e concentration. 

Recently I Whitehead et alj ( 1201 3h noted that the clusters which 
dissolve by the effects of stellar evolution lose their mass in a qual¬ 
itatively different way from those dominated by relaxation. 

The former, as they approach dissolution, lose the last frac¬ 
tion of their mass (which may be substantial) extremely rapidly, 
whereas the latter lose mass at a rate which is steady, and some¬ 
times even declining. Whitehead et al. also noted that the dividing 
line between the two modes of dissolution is quite sharp. For that 
reason it would not be surprising if a very small effect, such as 
the loss or retention of NS, were to place a cluster in one mode of 
dissolution or the other. _ 

Th e two kinds of behaviour described by IWhitehead et al.l 
( l2013h are plainly visi ble in several previous studies of star clus¬ 
ter evolution, such as iTakahashi & Portegies Zw^ (l2000l) . and 
they are visible in Fig. [T] of the present Letter, where all mod¬ 
els with kicks end their evolution by losing mass precipitately 
(except for the case N = 8k), whereas the others lose mass at a 
more moderate rate. We refer to these two cases as “jumping” and 
“skiing”, respectively. Fig. [T] also illustrates the point made by 
IChernoff & Weinberd ( 1990l) . i.e. the two modes of dissolution are 
characterised by the presence or absence of core collapse before 
dissolution. Indeed we see that the clusters with and without natal 
kicks (except for the case N = 8k) lie on either side of the divide 
between the two modes. 

In order to visualise the transition between skiing and jump¬ 
ing models, it has been particu larly instructive f or us to take on the 
point of view first suggested bv lWeinberd l ll993l) . and to explore the 
evolution of our models in the plane defined by the concentration 
(parameterized by c = log(rj/rc), where is the core radius) 
and the mass which remains bound to the system. In this repre¬ 
sentation, a system which experiences exclusively stellar evolution 
effects would gradually lose mass, while reducing its concentra- 


^ This unconventional terminology was coined by Simon Portegies Zwart 
in conversation with one of us (DCH) several years ago. It vividly conveys 
the differenc e between skiing down a gentle slope and jumping off a cliff. 
We note that I Whitehead et al] l2013h have conflated the two terms, with a 
different semantics. 


tion due to the progressive expansion, giving rise to a track movin g 
down in the plane and to the left (see Fig. [3 in IWeinberdll993h . 
The tracks qualitatively resemble those of some of our models, as 
shown in Fig.|4l These are four of the models with kicks, shown in 
red; and one of these (128kK) is also shown in Fig.[5l 

In Weinberg’s treatment, dealing with the slow evolution of 
spherical equilibrium models, the tracks end when equilibrium is 
no longer possible; the tracks end at points along a curve, which 
is shown as a dashed near-vertical curve in these figures. W-body 
models can cross this curve, but then lose mass on a dynamical time 
scale, explaining the jumping profile of the corresponding curves in 
Figs. [Handle Though its precise position may differ slightly when 
the simplifying assumptions of Weinberg’s models are relaxed, we 
refer to this curve as “Weinberg’s cliff”. 

In Weinberg’s models, mass-loss is driven by stellar evolution 
only, and his results should be applicable when this process dom¬ 
inates. When the effects of two-body relaxation are dominant, one 
of the natural consequences is the progressive increase of the cen¬ 
tral concentration, leading to core collapse. This results in a track 
oriented to the right-hand-side of the plane, behaviour which can 
be immediately recognised in the remaining models in Fig. |4] and 
It should not come as a surprise now that all long-lived, “skiing” 
models show signatures of core collapse, in contrast with short¬ 
lived, “jumping” models. 

These figures strongly suggest the existence of trajectories in 
which the two processes, stellar evolution and relaxation, are in a 
delicate balance overall, even though stellar evolution dominates 
early on and relaxation dominates thereafter. We have been partic¬ 
ularly fortunate to have included in our survey two models whose 
evolution almost perfectly delimits a “separatrix” between “skiing” 
and “jumping” systems (see the innermost pair of red lines in Fig.|4l 
which correspond to models 8kK and 16kK). Even more strikingly, 
the model 128kKe, despite the oscillations generated by the time- 
dependent tide, offers an excellent representation of the “separa¬ 
trix” (see the green line in Fig.^. 

Evidently, the models that we have studied lie close to the sep¬ 
aratrix dividing jumping models (which are dominated by stellar 
evolution, lose mass rapidly at the end of their lives, and do not 
reach core collapse) and skiing models (which are dominated by 
two-body relaxation, lose mass gently towards the end of their 
lives, and reach core collapse). If neutron s tars are given no natal 
kick, a s in model 128kN, or the models of iBaumgardt & Makind 
( l2003h . the trend to mass segregation and core collapse is accentu¬ 
ated, and the model moves across the separatrix into the domain of 
relaxation-dominated evolution. But we warn the reader against in¬ 
terpreting this as a general rule. Kicks were applied to both model 
8 kK and model 16kK (the innermost pair of red lines in Eig. 0, 
and they lie on opposite sides of the separatrix. The 8kK model, 
because of the low particle number and consequently smaller re¬ 
laxation time, is sufficiently dominated by relaxation to lie in the 
skiing regime. 

These considerations do not immediately explain, however, 
why the lifetime should be so different as a factor of nearly four. 
But the example of models 32kK and 8kN, which lose mass in al¬ 
most the same way until core collapse in the latter model (Eig[T]l, 
shows that the effects of skiing and jumping lead to different life¬ 
times. Though the difference is only a factor 1.13 in this case, it 
seems plausible that the effect could be much bigger if the event 
which determines the mode of dissolution occurred very early in 
the lifetime of a model, e.g. the ejection of neutron stars. Eurther- 
more, because our models lie so close to the separatrix between 
the two modes, it would not be surprising if very minor system- 
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c 

Figure 4. The plot shows the total mass remaining in the cluster as a func¬ 
tion of concentration, for the models depicted in Fig.^with (red lines) and 
without (blue lines) NS kicks. The black dots show the moment of core col¬ 
lapse, the black horizontal dashed line at M/Mq = 0.05 marks the formal 
dissolution condition, an d the black vert ical dashed line denotes “Wein¬ 
berg’s cliff” (see Fig. [^in lWeinberJ 1993) . 


atic differences in the initial conditions were to lead to significant 
systematic differences in the lifetime, as discussed in Sect. 3.2. 

While this Letter has focused on kicks by neutron stars, the les¬ 
son to learn is that apparently minor changes can have very large 
effects, especially for clusters close to the transition between differ¬ 
ent modes of dissolution. Other factors which should be taken into 
account include the presence and properties of primordial binaries, 
variations in the high-mass end of the IMF, and the degree of pri¬ 
mordial mass segregation, which influences both the importance of 
mass loss by stellar evolution and the role of remnants, not only NS 
but also stellar-mass black holes. The importance of these factors 
depends on the location of the dividing line between the two modes 
of dissolution that we have discussed, which can be assessed only 
by means of appropriate numerical experiments. 


4.2 Conclusions 

We have presented evidence, based on W-body simulations of the 
evolution of initially tidally filling King models with stellar evolu¬ 
tion, that the presence or absence of NS natal velocity kicks can 
play a crucial role in the long-term survival of model star clusters. 
In particular we show that some of the ba sic models in the land¬ 
mark study of lSaumgardt & Makind l l2003h are especially sensitive 
to this effect, which can change their lifetime by almost a factor of 
four. We explain this finding by showing that the models lie close to 
a dividing line between (i) models which are dominated by the ef¬ 
fects of mass-loss from stellar evolution, and whose evolution ends 
with a steepening rate of mass loss, and (ii) models whose dynam¬ 
ical evolution is dominated by two-body relaxation, which reach 
core collapse before dissolving, and do so with a gently decreasing 
rate of mass loss. 



Figure 5. As in Fig. |4] but presenting all models with N=128k. Note that 
model 128kKe (green line) spans the region occupied by the separatrix, dis¬ 
tinguishing “skiing” (128kN7, 128kNe, 128kN, 128kK7) and “jumping” 
models (128kK). 
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